scRNAseq data integration with Seurat
28 September 2022
0.1 A note on scRNAseq integration
Analysing single-cell datasets from differnet sources / conditions poses unique challenges. In particular, identifying cell populations that are present across multiple datasets can be problematic under standard workflows.
Seurat v4 includes a set of methods to match (or ‘align’) shared cell populations across datasets. These methods first identify cross-dataset pairs of cells that are in a matched biological state (‘anchors’), can be used both to correct for technical differences between datasets (i.e. batch effect correction), and to perform comparative scRNA-seq analysis of across experimental conditions.
THe scRNA-seq integration performed here is also described in Stuart*, Butler* et al, 2019. As an example of applicatin, see here: a comparative analysis of human immune cells (PBMC) in either a resting or interferon-stimulated state.
0.2 Integration goals
The following tutorial is designed to give you an overview of the kinds of comparative analyses on complex cell types that are possible using the Seurat integration procedure. Here, we address a few key goals:
- Create an ‘integrated’ data assay for downstream analysis
- Identify cell types that are present in both datasets
- Obtain cell type markers that are conserved in both control and stimulated cells
- Compare the datasets to find cell-type specific responses to stimulation
0.3 Setup the Seurat objects
The example data is in the SeuratData package.
0.3.1 Load the data
# install dataset
install.packages("https://seurat.nygenome.org/src/contrib/ifnb.SeuratData_3.0.0.tar.gz", repos = NULL,
type = "source")
# load dataset
ifnb <- LoadData("ifnb")0.3.2 Split into groups by condition (e.g. in vivo and ex vivo)
NB If you are using the scTransfrom method, follow methodology at bottom, and then perform visualisations etc
# split the dataset into a list of two seurat objects (stim and CTRL)
ifnb.list <- SplitObject(ifnb, split.by = "stim")
# normalise and identify variable features for each dataset independently
ifnb.list <- lapply(X = ifnb.list, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})
# select features that are repeatedly variable across datasets for integration
features <- SelectIntegrationFeatures(object.list = ifnb.list)0.4 Perform integration
0.4.1 Identify cellular anchors
We then identify anchors using the FindIntegrationAnchors() function, which takes a list of Seurat objects as input, and use these anchors to integrate the two datasets together with IntegrateData().
immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, anchor.features = features)0.4.2 Use anchors to generate an integrated data assay
# this command creates an 'integrated' data assay
immune.combined <- IntegrateData(anchorset = immune.anchors)0.5 Perform integrated analysis
Now we can run a single integrated analysis on all cells!
# specify that we will perform downstream analysis on the corrected data note that the
# original unmodified data still resides in the 'RNA' assay
DefaultAssay(immune.combined) <- "integrated"
# Run the standard workflow for visualization and clustering
immune.combined <- ScaleData(immune.combined, verbose = FALSE)
immune.combined <- RunPCA(immune.combined, npcs = 30, verbose = FALSE)
immune.combined <- RunUMAP(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindNeighbors(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindClusters(immune.combined, resolution = 0.5)0.6 Visualise results
Generate UMAP plots and label by meta data
(if this is in vivo and ex vivo would aim for little difference)
# Visualisation
p1 <- DimPlot(immune.combined, reduction = "umap", group.by = "stim")
p2 <- DimPlot(immune.combined, reduction = "umap", label = TRUE, repel = TRUE)
p1 + p2
To visualise the two conditions side-by-side, we can use the split.by argument to show each condition colored by cluster.
DimPlot(immune.combined, reduction = "umap", split.by = "stim")
0.7 Identify conserved cell type markers
To identify canonical cell type marker genes that are conserved across conditions, we provide the FindConservedMarkers() function. This function performs differential gene expression testing for each dataset/group and combines the p-values using meta-analysis methods from the MetaDE R package. For example, we can calculated the genes that are conserved markers irrespective of stimulation condition in cluster 6 (NK cells).
i.e. This method runs DE analysis on stimulated and stimulated separately, and finds which genes are high / low only in certain clusters. Then it sees which have significant p values in both stimulated and unstimulated version of that cluster. These are conserved markers
###Find conserved markers
# For performing differential expression after integration, we switch back to the original
# data
DefaultAssay(immune.combined) <- "RNA"
nk.markers <- FindConservedMarkers(immune.combined, ident.1 = 6, grouping.var = "stim", verbose = FALSE)
head(nk.markers)## CTRL_p_val CTRL_avg_log2FC CTRL_pct.1 CTRL_pct.2 CTRL_p_val_adj
## GNLY 0 2.259971 0.944 0.045 0
## FGFBP2 0 1.233024 0.505 0.020 0
## CLIC3 0 1.373076 0.597 0.024 0
## PRF1 0 1.016511 0.422 0.017 0
## CTSW 0 1.197380 0.531 0.029 0
## KLRD1 0 1.137006 0.495 0.019 0
## STIM_p_val STIM_avg_log2FC STIM_pct.1 STIM_pct.2 STIM_p_val_adj
## GNLY 0.000000e+00 2.2975368 0.957 0.059 0.000000e+00
## FGFBP2 4.224915e-162 0.7086055 0.259 0.015 5.937273e-158
## CLIC3 0.000000e+00 1.4152888 0.623 0.031 0.000000e+00
## PRF1 0.000000e+00 1.7526008 0.864 0.057 0.000000e+00
## CTSW 0.000000e+00 1.2863590 0.596 0.035 0.000000e+00
## KLRD1 0.000000e+00 1.2054121 0.552 0.027 0.000000e+00
## max_pval minimump_p_val
## GNLY 0.000000e+00 0
## FGFBP2 4.224915e-162 0
## CLIC3 0.000000e+00 0
## PRF1 0.000000e+00 0
## CTSW 0.000000e+00 0
## KLRD1 0.000000e+00 00.7.1 UMAP of conserved markers
We can explore these marker genes for each cluster and use them to annotate our clusters as specific cell types.
# plot expression of selected markers
FeaturePlot(immune.combined, features = c("CD3D", "SELL", "CREM", "CD8A", "GNLY", "CD79A", "FCGR3A",
"CCL2", "PPBP"), min.cutoff = "q9")
# use markers to assign cell type identities
immune.combined <- RenameIdents(immune.combined, `0` = "CD14 Mono", `1` = "CD4 Naive T", `2` = "CD4 Memory T",
`3` = "CD16 Mono", `4` = "B", `5` = "CD8 T", `6` = "NK", `7` = "T activated", `8` = "DC", `9` = "B Activated",
`10` = "Mk", `11` = "pDC", `12` = "Eryth", `13` = "Mono/Mk Doublets", `14` = "HSPC")
DimPlot(immune.combined, label = TRUE)
0.7.2 Dotplot of conserved markers
The DotPlot() function with the split.by parameter can be useful for viewing conserved cell type markers across conditions, showing both the expression level and the percentage of cells in a cluster expressing any given gene. Here we plot 2-3 strong marker genes for each of our 14 clusters.
# set factor levels of cellular identities
Idents(immune.combined) <- factor(Idents(immune.combined), levels = c("HSPC", "Mono/Mk Doublets",
"pDC", "Eryth", "Mk", "DC", "CD14 Mono", "CD16 Mono", "B Activated", "B", "CD8 T", "NK", "T activated",
"CD4 Naive T", "CD4 Memory T"))
# get a list of selected markers to plot
markers.to.plot <- c("CD3D", "CREM", "HSPH1", "SELL", "GIMAP5", "CACYBP", "GNLY", "NKG7", "CCL5",
"CD8A", "MS4A1", "CD79A", "MIR155HG", "NME1", "FCGR3A", "VMO1", "CCL2", "S100A9", "HLA-DQA1",
"GPR183", "PPBP", "GNG11", "HBA2", "HBB", "TSPAN13", "IL3RA", "IGJ", "PRSS57")
# plot
DotPlot(immune.combined, features = markers.to.plot, cols = c("blue", "red"), dot.scale = 8, split.by = "stim") +
RotatedAxis()
0.8 Identify differential expressed genes across conditions
0.8.1 Scatter plot to visualise DE genes in each cell type
Now that we’ve aligned the stimulated and control cells, we can start to do comparative analyses and look at the differences induced by stimulation. One way to look broadly at these changes is to plot the average expression of both the stimulated and control cells and look for genes that are visual outliers on a scatter plot. Here, we take the average expression of both the stimulated and control naive T cells and CD14 monocyte populations and generate the scatter plots, highlighting genes that exhibit dramatic responses to interferon stimulation.
e.g. Can compare in vivo and ex vivo lung cell types, and would expect gene expression to align. Check out the genes that don’t align, and see whether this matters or not.
theme_set(theme_cowplot())
# get average expression of T cells, in each condition
t.cells <- subset(immune.combined, idents = "CD4 Naive T")
Idents(t.cells) <- "stim"
avg.t.cells <- as.data.frame(log1p(AverageExpression(t.cells, verbose = FALSE)$RNA))
avg.t.cells$gene <- rownames(avg.t.cells)
# get average expression of CD14 monocytes, in each condition
cd14.mono <- subset(immune.combined, idents = "CD14 Mono")
Idents(cd14.mono) <- "stim"
avg.cd14.mono <- as.data.frame(log1p(AverageExpression(cd14.mono, verbose = FALSE)$RNA))
avg.cd14.mono$gene <- rownames(avg.cd14.mono)
# select genes to label
genes.to.label = c("ISG15", "LY6E", "IFI6", "ISG20", "MX1", "IFIT2", "IFIT1", "CXCL10", "CCL8")
# make plot
p1 <- ggplot(avg.t.cells, aes(CTRL, STIM)) + geom_point() + ggtitle("CD4 Naive T Cells")
p1 <- LabelPoints(plot = p1, points = genes.to.label, repel = TRUE)
p2 <- ggplot(avg.cd14.mono, aes(CTRL, STIM)) + geom_point() + ggtitle("CD14 Monocytes")
p2 <- LabelPoints(plot = p2, points = genes.to.label, repel = TRUE)
p1 + p2
As you can see, many of the same genes are upregulated in both of these cell types and likely represent a conserved interferon response pathway.
0.8.2 Find markers that are specifically expressed in cell types and conditions
Because we are confident in having identified common cell types across condition, we can ask what genes change in different conditions for cells of the same type. First, we create a column in the meta.data slot to hold both the cell type and stimulation information and switch the current ident to that column. Then we use FindMarkers() to find the genes that are different between stimulated and control B cells. Notice that many of the top genes that show up here are the same as the ones we plotted earlier as core interferon response genes. Additionally, genes like CXCL10 which we saw were specific to monocyte and B cell interferon response show up as highly significant in this list as well.
# make a new column with celltype and condition
immune.combined$celltype.stim <- paste(Idents(immune.combined), immune.combined$stim, sep = "_")
immune.combined$celltype <- Idents(immune.combined)
# make this an identity
Idents(immune.combined) <- "celltype.stim"
# find markers unqieu to that identity
b.interferon.response <- FindMarkers(immune.combined, ident.1 = "B_STIM", ident.2 = "B_CTRL", verbose = FALSE)
head(b.interferon.response, n = 15)## p_val avg_log2FC pct.1 pct.2 p_val_adj
## ISG15 8.657899e-156 1.7770568 0.998 0.240 1.216695e-151
## IFIT3 3.536522e-151 1.9274217 0.964 0.052 4.969874e-147
## IFI6 1.204612e-149 1.8534970 0.966 0.080 1.692841e-145
## ISG20 9.370954e-147 0.9605165 1.000 0.672 1.316900e-142
## IFIT1 8.181640e-138 1.8102156 0.912 0.032 1.149766e-133
## MX1 1.445540e-121 1.5587694 0.907 0.115 2.031417e-117
## LY6E 2.944234e-117 1.4617133 0.894 0.152 4.137531e-113
## TNFSF10 2.273307e-110 1.6182250 0.792 0.025 3.194678e-106
## IFIT2 1.676837e-106 1.5821146 0.786 0.035 2.356459e-102
## PLSCR1 3.279290e-94 1.3253653 0.797 0.117 4.608387e-90
## IRF7 1.475385e-92 1.2257943 0.838 0.190 2.073358e-88
## CXCL10 1.350777e-82 1.7099596 0.641 0.010 1.898247e-78
## UBE2L6 2.783283e-81 1.0179545 0.851 0.300 3.911348e-77
## PSMB9 2.638374e-76 0.7067650 0.941 0.573 3.707707e-72
## EPSTI1 3.387863e-76 1.2033915 0.715 0.105 4.760964e-720.8.3 Feature plots and Violin plots
Its useful to visualise changes in gene expression with the split.by option of the FeaturePlot() or VlnPlot() function.
Feature plots
Display FeaturePlots of a list of given genes, split by a grouping variable (e.g. in vivo versus ex vivo).
In the example here:
CD3D and GNLY are canonical cell type markers (for T cells and NK/CD8 T cells) that are virtually unaffected by interferon stimulation and display similar gene expression patterns in the control and stimulated group.
IFI6 and ISG15, are core interferon response genes and are upregulated in all cell types.
CD14 and CXCL10 are genes that show a cell type specific interferon response. CD14 expression decreases after stimulation in CD14 monocytes, CXCL10 shows a distinct upregulation in monocytes and B cells after interferon stimulation but not in other cell types.
FeaturePlot(immune.combined, features = c("CD3D", "GNLY", "IFI6"), split.by = "stim", max.cutoff = 3,
cols = c(lightergrey, colors[1])) + t2
# NB Hack function if i want to change aesthetics
plots <- VlnPlot(immune.combined, features = c("LYZ", "ISG15", "CXCL10"), split.by = "stim", group.by = "celltype",
pt.size = 0, combine = FALSE)
wrap_plots(plots = plots, ncol = 1)
0.9 Save dataset
1 Performing integration on datasets normalised with SCTransform
In Hafemeister and Satija, 2019, we introduced an improved method for the normalisation of scRNA-seq, based on regularised negative binomial regression. The method is named ‘sctransform’, and avoids some of the pitfalls of standard normalization workflows, including the addition of a pseudocount, and log-transformation. You can read more about sctransform in the manuscript or our SCTransform vignette.
Below, we demonstrate how to modify the Seurat integration workflow for datasets that have been normalised with the sctransform workflow. The commands are largely similar, with a few key differences:
- Normalise datasets individually by
SCTransform(), instead ofNormalizeData()prior to integration - As discussed further in our SCTransform vignette, we typically use 3,000 or more features for analysis downstream of sctransform.
- Run the
PrepSCTIntegration()function prior to identifying anchors - When running
FindIntegrationAnchors(), andIntegrateData(), set thenormalization.methodparameter to the valueSCT. - When running sctransform-based workflows, including integration, do not run the
ScaleData()function
LoadData("ifnb")
ifnb.list <- SplitObject(ifnb, split.by = "stim")
ifnb.list <- lapply(X = ifnb.list, FUN = SCTransform)
features <- SelectIntegrationFeatures(object.list = ifnb.list, nfeatures = 3000)
ifnb.list <- PrepSCTIntegration(object.list = ifnb.list, anchor.features = features)immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, normalization.method = "SCT",
anchor.features = features)
immune.combined.sct <- IntegrateData(anchorset = immune.anchors, normalization.method = "SCT")immune.combined.sct <- RunPCA(immune.combined.sct, verbose = FALSE)
immune.combined.sct <- RunUMAP(immune.combined.sct, reduction = "pca", dims = 1:30)p1 <- DimPlot(immune.combined.sct, reduction = "umap", group.by = "stim")
p2 <- DimPlot(immune.combined.sct, reduction = "umap", group.by = "seurat_annotations", label = TRUE,
repel = TRUE)
p1 + p2
Now that the datasets have been integrated, you can follow the previous steps in this vignette identify cell types and cell type-specific responses.
Session Info
sessionInfo()## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur ... 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] cowplot_1.1.1 metap_1.8 multtest_2.52.0
## [4] Biobase_2.56.0 BiocGenerics_0.42.0 ggplot2_3.3.6
## [7] patchwork_1.1.2 ifnb.SeuratData_3.0.0 SeuratData_0.2.2
## [10] sp_1.5-0 SeuratObject_4.1.2 Seurat_4.2.0
## [13] BiocStyle_2.24.0
##
## loaded via a namespace (and not attached):
## [1] sn_2.1.0 plyr_1.8.7 igraph_1.3.5
## [4] lazyeval_0.2.2 splines_4.2.1 listenv_0.8.0
## [7] scattermore_0.8 qqconf_1.3.0 TH.data_1.1-1
## [10] digest_0.6.29 htmltools_0.5.3 magick_2.7.3
## [13] fansi_1.0.3 magrittr_2.0.3 tensor_1.5
## [16] cluster_2.1.4 ROCR_1.0-11 globals_0.16.1
## [19] matrixStats_0.62.0 sandwich_3.0-2 spatstat.sparse_2.1-1
## [22] colorspace_2.0-3 rappdirs_0.3.3 ggrepel_0.9.1
## [25] rbibutils_2.2.9 xfun_0.33 dplyr_1.0.10
## [28] crayon_1.5.1 jsonlite_1.8.0 progressr_0.11.0
## [31] spatstat.data_2.2-0 survival_3.4-0 zoo_1.8-11
## [34] glue_1.6.2 polyclip_1.10-0 gtable_0.3.1
## [37] leiden_0.4.3 future.apply_1.9.1 abind_1.4-5
## [40] scales_1.2.1 mvtnorm_1.1-3 DBI_1.1.3
## [43] spatstat.random_2.2-0 miniUI_0.1.1.1 Rcpp_1.0.9
## [46] plotrix_3.8-2 viridisLite_0.4.1 xtable_1.8-4
## [49] reticulate_1.26 spatstat.core_2.4-4 stats4_4.2.1
## [52] htmlwidgets_1.5.4 httr_1.4.4 RColorBrewer_1.1-3
## [55] TFisher_0.2.0 ellipsis_0.3.2 ica_1.0-3
## [58] farver_2.1.1 pkgconfig_2.0.3 sass_0.4.2
## [61] uwot_0.1.14 deldir_1.0-6 utf8_1.2.2
## [64] labeling_0.4.2 tidyselect_1.1.2 rlang_1.0.6
## [67] reshape2_1.4.4 later_1.3.0 munsell_0.5.0
## [70] tools_4.2.1 cachem_1.0.6 cli_3.4.1
## [73] generics_0.1.3 mathjaxr_1.6-0 ggridges_0.5.4
## [76] evaluate_0.16 stringr_1.4.1 fastmap_1.1.0
## [79] yaml_2.3.5 goftest_1.2-3 knitr_1.40
## [82] fitdistrplus_1.1-8 purrr_0.3.4 RANN_2.6.1
## [85] pbapply_1.5-0 future_1.28.0 nlme_3.1-159
## [88] mime_0.12 formatR_1.12 compiler_4.2.1
## [91] rstudioapi_0.14 plotly_4.10.0 png_0.1-7
## [94] spatstat.utils_2.3-1 tibble_3.1.8 bslib_0.4.0
## [97] stringi_1.7.8 highr_0.9 rgeos_0.5-9
## [100] lattice_0.20-45 Matrix_1.5-1 vctrs_0.4.1
## [103] mutoss_0.1-12 pillar_1.8.1 lifecycle_1.0.2
## [106] BiocManager_1.30.18 spatstat.geom_2.4-0 Rdpack_2.4
## [109] lmtest_0.9-40 jquerylib_0.1.4 RcppAnnoy_0.0.19
## [112] data.table_1.14.2 irlba_2.3.5 httpuv_1.6.6
## [115] R6_2.5.1 bookdown_0.29 promises_1.2.0.1
## [118] KernSmooth_2.23-20 gridExtra_2.3 parallelly_1.32.1
## [121] codetools_0.2-18 MASS_7.3-58.1 assertthat_0.2.1
## [124] withr_2.5.0 mnormt_2.1.0 sctransform_0.3.5
## [127] multcomp_1.4-20 mgcv_1.8-40 parallel_4.2.1
## [130] grid_4.2.1 rpart_4.1.16 tidyr_1.2.1
## [133] rmarkdown_2.16 Rtsne_0.16 numDeriv_2016.8-1.1
## [136] shiny_1.7.2